
Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 


1. Thesis and Dissertation Collection, all items 


1996-12 

Ocean wave data analysis using Hilbert 
transform techniques 


Navarro, Moises M. 

Monterey, California. Naval Postgraduate School 


http://hdl.handle.net/10945/32022 


Downloaded from NPS Archive: Calhoun 



DUDLEY 

KNOX 

LIBRARY 


http://www.nps.edu/ljbrary 


CsMwun is the Naval Postgraduate School's public access distal repository for 
research mate riels and tnstitutiooal pubfications created by the NPS community. 
Cathoufii is named for Professor of Mathematics Guy K. CatHiuo, NPS's first 
appoinited — and publi^d — scholar^ author. 

Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 
MontereVr California USA 93943 



NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CALIFORNIA 



THESIS 



OCEAN WAVE DATA ANALYSIS USING HILBERT 
TRANSFORM TECHNIQUES 

by 

Moises M. Navarro 
December, 1996 

Thesis Advisor: Andres Larraza 

Co-Advisor: Roberto Cristi 


Approved for public release; distribution is unlimited. 

199T0919 033 




REPORT DOCUMENTATION PAGE 

Form Approved 0MB No. 0704-0188 

, f- 

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instruction, searching 
existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this 
burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, 
Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and 
Budaet. Paperwork Reduction Project (0704-0188) Washington DC 20503. 

1. AGENCY USE ONLY 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 

December, 1996. Master’s Thesis 

4. TITLE AND SUBTITLE 

OCEAN WAVE DATA ANALYSIS USING HILBERT TRANSFORM 
TECHNIQUES 

5. FUNDING NUMBERS 

6. AUTHOR Moises Navarro Lopez 

Lieutenant Commander, Venezuelan Navy 

7. PERFORMING ORGANIZATION NAME AND ADDRESS 

Naval Postgraduate School 

Monterey CA 93943-5000 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

Naval Postgraduate School 

Monterey CA 93943-5000 

10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

11 . SUPPLEMENTARY NOTES The views expressed in this thesis are those of the author and do not reflect the 
official policy or position of the Department of Defense or the U.S. Government. 

12a. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution is unlimited. 

12b. DISTRIBUTION CODE 

13. ABSTRACT 

A novel technique to determine the phase velocity of long-wavelength shoaling waves is 
investigated. Operationally, the technique consists of three steps. First, using the Hilbert transform 
of a time series, the phase of the analytic signal is determined. Second, the correlations of the 
phases of analytic signals between two points in space are calculated and an average time of travel 
of the wave fronts is obtained. Third, if directional spectra are available or can be determined from 
time series of large array of buoys, the angular information can be used to determine the true time of 
travel. The phase velocity is obtained by dividing the distance between buoys by the correlation time. 
Using the Hilbert transform approach, there is no explicit assumption of the relation between 
frequency and wavenumber of waves in the wave field, indicating that it may be applicable to 
arbitrary wave fields, both linear and nonlinear. Limitations of the approach are discussed. 

14. SUBJECT TERMS 

Hilbert Transform, Sea waves. Cross-correlations 

15. NUMBER OF PAGES 

52 

16. PRICE CODE 

17. SECURITY CLASSIFICA¬ 
TION OF REPORT 

Unclassified 

18, SECURITY CLASSIFI¬ 
CATION OF THIS PAGE 

Unclassifieid 

19. SECURITY CLASSIFICA¬ 
TION OF ABSTRACT 

Unclassified 

20. LIMITATIONOF 
ABSTRACT 

UL 


NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89) 


Prescribed by ANSI Std. 239-18 298-102 


1 
















11 




Approved for public release; distribution is unlimited. 


OCEAN WAVE DATA ANALYSIS USING HILBERT TRANSFORM 

TECHNIQUES 

Moises M. Navarro 

Lieutenant Commander, Venezuelan Navy 
B.S., Venezuelan Naval School, 1982 

Submitted in partial fulfillment 
of the requirements for the degree of 

MASTER OF SCIENCE IN APPLIED PHYSICS 

from the 


Author: 


Approved by: 




iii 






IV 




ABSTRACT 


A novel technique to determine the phase velocity of long-wavelength 
shoaling waves is investigated. Operationally, the technique consists of three 
steps. First, using the Hilbert transform of a time series, the phase of the 
analytic signal is determined. Second, the correlations of the phases of analytic 
signals between two points in space are calculated and an average time of travel 
of the wave fronts is obtained. Third, if directional spectra are available or can 
be determined from time series of large array of buoys, the angular information 
can be used to determine the true time of travel. The phase velocity is obtained 
by dividing the distance between buoys by the correlation time. Using the Hilbert 
transform approach, there is no explicit assumption of the relation between 
frequency and wavenumber of waves in the wave field, indicating that it may be 
applicable to arbitrary wave fields, both linear and nonlinear. Limitations of the 
approach are discussed. 





TABLE OF CONTENTS 


I. INTRODUCTION.1 

II. SURFACE WAVES IN THE HAMLET’S COVE I EXPERIMENT.5 

A. SURFACE WAVES IN WATER.5 

B. HAMLET’S COVE TOPOGRAPHY.7 

III. THE HILBERT TRANSFORM.11 

IV. DATA ANALYSIS....17 

A. THECHEBYSHEV TYPE IIFILTER.17 

B. THE XCORRELATION FUNCTION.19 

C. PROCEDURE AND RESULTS.19 

V. DISCUSSION OF RESULTS AND RECOMMENDATIONS FOR FUTURE 

WORK . 27 

APPENDIX .. 29 

LIST OF REFERENCES.37 

INITIAL DISTRIBUTION LIST.41 


vii 


















ACKNOWLEDGMENTS 


First and foremost I would like thanks God for helping me to overcome 
one more stage in my life that was longer and difficult than I though. 
Particularly to my friend, lover and spouse Carmen Luisa and our children, 
Moises, David and Daniel; they with my mother Herlinda enduring me 
throughout this long course. Their love, support and encouragement helped to 
study at the Naval Postgraduate School, and they were always my inspiration. 
All of you make this experience to have a happy end, thank you guys. 

My sincere thanks to Prof. Andres Larraza who gives me the opportunity 
to work with him in this technique and helped me prepare and Implement this 
Thesis. His guidance, patience, and insight enabled me to learn a tremendous 
amount in a very short period of time. Thanks Andres and God bless you. 

I also want to express in these lines my grateful to Prof. Robert CristI who 
without any inconvenient and with all his enthusiasm help sincerely In the 
accomplishment of this thesis. 

Finally I would like to thanks in the Combat System Curriculum Comdr. 
Michael Witt and Prof. James Sander. They helped and guided me to work out 
and finding the way to make this master degree a successful experience. 





I. INTRODUCTION 


As ocean waves shoal, the wave field evolves substantially from Its deep 
water state. While the geometry of the dynamic currents and of the bottom 
topography determine the direction and to some extent the energy of the spectral 
components, nonlinearities in the wave field lead to harmonic generation and to 
spectrum broadening. Furthermore, for a nonlinear wave field, there is an 
apparent negative refraction of the spectral peak (Abreu et al. 1992). Thus, for a 
known sloping bottom, the refraction away from the beach normal is a direct 
consequence of nonlinear energy transfer. On the other hand, linear theory 
does not predict most of these changes, and for a given shoaling spectrum, 
linear theory might lead to an incorrect determination of the beach normal. 

Directional spectrum has become commonly used to describe nonlinear 
random ocean waves because it provides ways to estimate the significant wave 
heights in a breaker and also to determine the predominant direction for net 
sediment transport. Considerable progress has been made over the last two 
decades to study the evolution of the wave spectrum due to refraction, 
diffraction, and nonlinear interaction under conditions of mild slope (Lui et al. 
1985, Freilich and Guza 1984, Abreu et al. 1992). However, field conditions 
often violate the basic assumption that the wavelength be smaller than the 
characteristic length scale of bottom variations. Nevertheless, it can be shown 
(Larraza, 1994) that using a ray optics approximation, the determination of the 
kinematic wave parameters (wavelength, direction of propagation) for a given 
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geometry is within 5% of its true value. On the other hand, the determination of 
the dynamic values (energy, momentum) exceeds at least 15% the true value. 
Thus by using ray theory, sediment transport and other relevant effects caused 
by the wave’s momentum will be ill estimated. 

For an arbitrary bottom topography, wave energy and momentum can be 
calculated within the physical optics approximation. The passage from 
geometrical to physical optics can in principle be accomplished by use of a path 
integral formulation approach. Within this formulation, the basis states are 
characterized by the "rays", and the physical optics results from a sum over all 
the possible paths (and not only the path determined by Fermat’s principle). The 
path integral formulation approach has not been considered previously either for 
gravity waves over a changing topography or for shoaling waves. This situation 
is particularly surprising because the problem can be stated from the onset by 
using the known results of ray optics, either linear or nonlinear, and performing a 
sum over all possible ray paths of the basis states. Even though there are 
standard analytical techniques to deal with certain special cases, the most 
general ones can only be solved numerically. 

Because the sea floor is known to be a controlling factor in low-frequency 
shallow water wave propagation (Akal, 1980, Akal and Jensen 1983), phase 
velocity measurements can in principle be used to determine bottom topography, 
which can in turn be used for true ground determination in remote sensing 
applications. The inverse problem, that is, the determination of the topography 
given the evolution of the spectrum along a path requires an accurate technique 
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to determine celerity and direction of the waves in a given wave field data. 
Common techniques of Fourier transform are usually inadequate. In the Fourier 
transform, a real space-time signal is converted to a complex frequency- 
wavevector signal. Thus, for applications to ocean wave field data a dispersion 
relation between frequency and wavenumber has to be assumed in order to 
determine the phase velocity of the spectral components. The problem can be 
complicated by nonlinearities in the wave field. 

On the other hand, the Hilbert transform of a real valued time-space signal 
is another real value time-space signal for which the phase of a signal can be 
calculated in terms of the signal Itself and its Hilbert transform. Thus by using 
the time series of a signal measured at different points, one can determine a 
relative phase shift and extract phase velocity measurements without reference 
to a dispersion relation. 

In this thesis we present an analysis of ocean surface waves data 
obtained by The Naval Research Laboratory Stennis Center, in support of 
Hamlet’s Cove I (Smith, 1994). An Immediate goal is to extract phase velocity 
information from time series at three buoys. The data analysis technique used 
throughout this thesis is based on the Hilbert transform. This technique allow us 
to get local properties, local energy and local frequency, of the signal instead of 
the global properties that others techniques, like the Fourier transform, would 
give (Long, 1995). 

In Chapter II we give a brief introduction to linear water waves and point 
out their dispersive character. We also present a description of the local 
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topography for the wave data that we analyze and emphasize the need for an 
inverse problem in arbitrary topographies. Chapter III gives a detailed account 
of the basic ideas involved in Hilbert transforms. In Chapter IV we present the 
results and the procedures used for the analysis of the data. The data was 
analyzed using the Fast Fourier Transform (FFT), the Hilbert transform, and the 
cross-correlation functions of MATLAB. In Chapter V we discuss the results and 
provide some recommendations for future work. 
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II. SURFACE WAVES IN THE HAMLET’S COVE I EXPERIMENT 


In the first part of this chapter, we present a brief introduction to surface 
gravity waves over uniform depth and emphasize their dispersive behavior. In 
contrast, the data that we analyze throughout this thesis is over nonuniform 
terrain. The second part of this chapter gives a detailed description of the local 
topography for the Hamlet’s Cove I experiment. 

A. SURFACE WAVES IN WATER 

When a group of waves moves across the surface of water, each 
particular wavecrest travels faster than the group as a whole, and eventually 
passes through it. Thus new crest continually are being created at the back of 
the group while old crests are disappearing at the front. 

This behavior does occur because water waves are dispersive. We can 
say that this implies that the different Fourier components that make the general 
disturbance have phase velocities that depend on their wavelength. For surface 
waves, if the surface elevation of water with uniform depth h is described by the 
wavetrain 

Tj = Acos(kx- cot) , (1.1) 

the wave speed is given by 
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c = CO / k 


^ tanh(kh) 

Iv. 


\l/2 

J 


( 1 . 2 ) 


so that waves of longer wavelength, X = 2K/k, travel faster. In Eq (1.1) k is the 
wavenumber and co is the frequency of the wave. 

While each individual wavecrest travels with speed c, the velocity of travel 
of the group as a whole is the group velocity 

Cg = doo/dk . (1.3) 

The phase velocity in water with uniform depth h (Acheson, 1995) can not 
exceed Vgh. For kh large, i.e. h »X, tanh(kh) = 1 , and 


= g / k , (1.4) 

corresponding to the case of infinite depth. A good approximation for deep water 
waves is for depths h greater than about ^X. In shallow water, i.e. h«XI 2k, 

tanh(kh) = kh , and equation (1.2) becomes 



(1.5) 
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This equation is important, because as we can see c is independent of k for 
shallow water. As we will see, the data analyzed in this thesis is remarkably 
nondispersive. 

B. HAMLETS COVE TOPOGRAPHY 

The data analyzed in this thesis is from three bottom-mounted pressure 
gauge/ current meters for the measurements correspond to directional wave 
spectra, currents, and tidal elevation. The gauges were deployed in depths of 
10’, 15’, and 25’ installed in support of Hamlet Cove I experiment. 

The three bottom-mounted pressure gauges were placed on either side of 
the sand bar as shown in Figures 2.1 and 2.2 

86<>48’10” 86“48’20” 

30<>23’10” 

30“23’00” 

30°22’50” 

Figure 2.1. Coordinates for the Positions of the gauges and their depths in feet. 



• #3 (Depth 10’) 

« #2(Depth 15’) 


#1 (Depth 25’r 
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Figure 2.2. Distance from shore (m) of the three buoys and reference to the sand bar. The sand 
bar profile, measured relative to the surface, is a t a depth of 1.3 m 150 m off shore. 

The buoys used collected data from the wave as pressure, x-component 
of the velocity and the y-component as well. The sampling was continuous with 
a sampling frequency of 2 Hz. 

Because the pressure gauge responds to the rise and fall to the free 
surface, the pressure at depth z attenuates by a factor proportional to e ' 

where k is 2TdX and X is the wavelength of the Fourier component of the surface 
disturbance. At the 15’ and 25’ depths, (gauges 1 and 2) frequencies above 0.3 
Hz were below the noise level of the instruments and, therefore, were ignored in 
the analysis. The high frequency cut-off for gauge number 3 is 0.5 Hz. 

Table 2-1 shows the data collection periods for each instrument. We analyzed 
the data collected by the three buoys in 16 July 1994. 
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Gauge 

1 

2 

3 

DATES 

14 JUL-20JUL 

14 JUL-20 JUL 

14 JUL- 20 JUL 

20 JUL-29 JUL 

20 JUL- 29 JUL 

20 JUL- 29 JUL 

29 JUL-1 AUG 

29 JUL- 3 AUG 

29 JUL- 3 AUG 

3 AUG-10 AUG 

3 AUG-10 AUG 

3 AUG-9 AUG 

12 AUG-17 AUG 

10 AUG-17 AUG 

FAILED 

18 AUG-24 AUG 

18 AUG-24 AUG 

17AUG-24 AUG 

24 AUG-4 SEP 

FAILED 

FAILED 


Table 2-1 Dates for the data collected at the three buoys. The data analyzed in this thesis 
corresponds to data collected on 16 July 1994. 


From the bathymetry survey shown in Figure 2.2, we can see that 
refraction and diffraction effects from waves 60 m long and longer may be used 
to characterize the bottom topography. This long waves, due mainly to swell, are 
essentially nondispersive and already violate the conditions for ray optics to 
apply. Because shorter dispersive waves are mainly due to local winds, they 
would not give accurate topography estimates. Thus, phase velocity estimates 
from the long waves are the first step to determine an unknown bathymetry. 
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III. THE HILBERT TRANSFORM 


In this chapter we review some of the basic concepts of the Hilbert 
transform, following the review articie by Bendat (1989). In order to elucidate the 
concept and possible applications of the Hilbert transform, we present three 
equivalent definitions, namely as a convolution integral, as a 7u/2 phase shifter, 

and as an imaginary part of an analytic signal. 

The Hilbert Transform of a real valued signal is a complex signal called 
the analytic signal. The real part of the analytic signal is the original real-valued 
time signal while the imaginary part is a copy of the original signal with each of its 
Fourier components shifted in phase by 90 (Bendat, 1989). The transformed 
signal retains the same amplitude and frequency information contained in the 
original signal with the added phase information dependent on the phase of the 
original data. 

For any real time series, C(t), the Hilbert transform ^ (t) is 



in which P implies the principal value. The analytic continuation of the Hilbert 
transform is thus the Cauchy integral which corresponds to translations of an 
analytic function in the complex domain of integration. 

The analytic signal Z(t) is the real based time series defined by 
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z (t) = m + j 4(t) 


(3.2) 


which can also be written as 

2(t) = A(t) e , (3.3) 

where A(t) is called the envelope signal or the amplitude and i&(t) is called the 
instantaneous phase signal, defined by: 


A(t) = [?"(t) + |^(t)]“ 


(3.4) 


and 


«(t) = tan-'[«t) / «t)], 


(3.5) 


respectively. The instantaneous frequency co is given by 


CD 


d^(t) 

dt 


(3.6) 


and the local energy is one half of the square value of the amplitude. Eqs. (3.4), 
(3.5), and (3.6) represents local properties in time. 

The properties of the Hilbert Transform can be best understood in terms 
of three equivalent definitions. 

1) Definition as Convolution Integral 

The Hilbert Transform of a real valued function ^(t) extending from - ~ < t 
< + ~ is a real valued function as described in (3.1): 
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(3.7) 




thus ^ (t) is the convolution integral of C(f) with (1 / tc t), written as 


/ H ^ 


ut) = C (t) 


V7t ty 


(3.8) 


2) Definition as (tc / 2) Phase Shift System 


If F'{ f) is the Fourier transform of ^ (t), namely 


TOO 

F'(<) = JS (t)e-'^”"dt, 


(3.9) 


then, from the convolution property of the Hilbert transform (Bendat, 1989) and 
equation (3.7), it follows that F^( f) is the Fourier transform of C(t), multiplied by 

the Fourier transform of (1 / ic t) where 


7 e 

—dt = -jsgnf 


-j for f > 0 
j for f < 0 


(3.10) 


Hence, equation (3.6) is equivalent to the passage of ^ (t) through a system 
defined by (-j sgnf) to yield: 

F'(f) = (-isgnf)F"(f) , (3.11) 
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where F^^(f) is the Fourier transform of ^(t). The complex-valued quantity F^( f) is 
the Hilbert transform of the complex-valued quantity F^^ (f) as defined by: 


F'(f) = HF'/(f)]= (-jsgnf)F'^(f) 


/// 


The Fourier transform (-j sgnf) can be represented by 


(3.12) 


B (f) = - j sgnf = 


[e forf > 0 I 
e ^ ^ for f < 0 


(3.13) 


or using the complex polar notation 


r I -j<l> (f) 

B(f) = |B(f)| e ■> , 


(3.14) 


Hence, B(f) is a (tc / 2 ) phase shift system where 


|B(f)|= 1 for all f, 


(3.15) 


and 


<P b (f) = 


[ Tt / 2 for f > 0 1 
I-7C/2 for f<0r 


(3.16) 


Thus if one writes: 


(f) = \ f " (f)| , 


(3.17) 
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it follows that 


F’(f) = |F’(f)| = |F''(f)| , (3.18) 

Thus the Hilbert transform consist of passing C(t) through a system which leaves 
the magnitude of f" (f) unchanged, but changes the phase from <}) x (f) to (j) x (f) + 
(|) b (f). Thus 

9 x(f) ^ 9 X (0 + (’t / 2) for f > 0 , (3.19a) 

9x(f) ^ ‘Px(f) - (7t / 2) for f < 0 (3.19b) 

in other words, the Hilbert transform shifts by Ji / 2 for positive frequencies and 
by -n/2 for negative frequencies. 

3) Definition as Imaginary part of Analytic Signal 

A third useful way to understand and to compute the Hilbert transform is 
to introduce the analytic signal Z(t) which allows us to work with the amplitude 
A(t), and the phase ^ (t) equations (3.4 and 3.5). Time derivatives of the phase 

yield the instantaneous frequency co (3.6). 

Examples of Hilbert transforms and the corresponding Fourier transform 
of some simple functions are shown in Figure 3.1. The convolution and phase 
shifts due to the Hilbert transform are apparent. Also apparent is the fact that 
the analytic signal preserves local amplitude. 
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Figure 3.1. Some examples of Hilbert Transform. 
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IV. DATA ANALYSIS 


In this chapter we present the procedure and the results of the analysis of 
the Hamlet’s Cove I surface wave data. The main goal is to determine if there is 
any phase correlation between data at the three buoys after low pass filtering 
and applying the Hilbert transform to the filtered data. Phase correlations can 
provide valuable information because they can tell us if the signal recorded from 
the three buoys corresponded to the same wave and if so, they can yield directly 
the phase velocity of the wave. Phase velocity is obtained by interpolating the 
correlation time between the buoys for which the location is known. 

The data analyzed is of pressure measurements from a pressure gauge 
and two orthogonal components of velocity from a current meter. Because 
pressure gauge measurements are more sensitive to long waves, we will 
consider the low pass filtered output from for pressure gauge of measurements 
on 16 July 1994. 

A. THE CHEBYSHEV TYPE II FILTER 

An optimal filter needs to maintain amplitude information with no phase 
distortion in the pass band. In our case, the latter requirement is particularly 
important because from our approach, phase velocity values can only be 
determined from phase correlations. Because the filtering process should 
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minimize distortion of all small amplitude fluctuations, the stop band has to be at 
least 40 dB below the pass band to eliminate possible leakage problems and it 
has to have a minimum transition width to ensure good frequency separation. 
Most of these requirements can be met by with a Chebyshev type II filter. Higher 
order filters could introduce numerical errors, so the filter has to meet the 
requirements with the lowest order possible. 

As shown in the Appendix (mfile), we used Matlab Signal Processing 
Toolbox (Mathworks, 1994). Of the two techniques, Analog Prototyping and 
Direct Design, we chose Analog Prototyping because the Direct Design 
technique has very stringent numerical accuracy constraints. 

The Analog Prototyping technique allows the construction of digital 
equivalents of certain classical analog filters, especially Butterworth, Bessel, 
Chebyshev, and Elliptic. Among them and for ours purpose the Chebyshev Type 
II filter has the narrower transition and the flatter pass band response needed for 
the analysis of the data. 

Analog filters are HR (infinite-duration impulse) type filters, which 
necessarily introduce phase distortion. In order to overcome this problem the 
data was filtered both directions, forward and reverse by using the Matlab 
function filtfilt (Mathworks, 1994). With this approach, we can get precise zero- 
phase distortion and doubling of the filter order. Also, filtfilt minimizes startup 
and ending transients by adjusting initial conditions. 
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B. THE XCORRELATION FUNCTION 

Matlab’s function xcorr estimates the cross-correlation sequence of 
random processes. For the purpose of determining phase velocity values from 
apparently random time series, the cross-correlation function is important 
because it can provide a measure of the similarity between two signals. 

The xcorr function can calculate the cross-correlation between two arrays 
of vectors of the same length, the autocorrelation for a vector, and the cross¬ 
correlation for a matrix. In all these forms of xcorr, the zero lag of the output is in 
the middle of the sequence. In order to normalized the sequence we used the 
option ‘coeff so the autocorrelations at zero lag are identically 1.0. 

C. PROCEDURE AND RESULTS 

Figure 4.1 shows the Fourier transform of the pressure time series for the 
three buoys. From the shape of the spectrum, we can conclude that waves due 
to the local wind may correspond to frequencies above 0.1 Hz. If phase velocity 
values can be used to determine local bathymetry, waves generated by the local 
wind are unwelcome noise for the purpose of our analysis. Instead, the longer 
waves in the swell may characterize the bathymetry more accurately because of 
their direct response to the bottom topography manifested by refraction and 
diffraction effects. Swell may correspond to frequencies below 0.08 Hz, which 
we chose as the cutoff frequency for the low pass band of the filter. 
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Frequency (Hz) 




Figure 4.1. FFT of the pressure data for the three buoys. Frequencies above 0.1 Hz may 

correspond to waves generated by the local wind as evidenced by the increase of 
spectral energy and range. 

With this information in mind, we filtered the whole time series using a 
Chebyshev Type II filter with a cut off frequency of 0.08 Hz. Because of the 
Chebyshev Type II filter introduces phase distortion, the wave data was filtered 
in both the forward and reverse directions using the MATLAB function filtfilt, thus 


20 






obtaining precise zero-phase distortion and doubling of the filter order. As 
mentioned before, filtfilt minimizes startup and ending transients by adjusting 
initial conditions. 

We applied the Hilbert transform to the filtered time series. Figure 4.2 
shows the Matlab Power Spectrum function (PSD) of unfiltered time series and 
the PSD of the Hilbert transform of the data from the first buoy after filtering. 
Apparent from the figure, is that the power spectrum of the unfiltered data does 
not exhibit very good frequency-energy decomposition. On the other hand, the 
power spectrum of the Hilbert transform of the filtered time series shows good 
frequency-energy decomposition and enhancement of the more energetic waves 
of the filtered time series. 

The data from each buoy has 172,800 points corresponding to two data 
point per second in a full day. Because this large array is numerically very 
demanding, we chose to work with a window of two hours, or 10% of the whole 
data. We applied Matlab’s Hilbert transform (^(t)) function to the filtered data in 
this time window. Along with the real time series (Y1), MATLAB generates the 
analytical signal 

Z(t) = Y1(t) + i^ (t) , (4.1) 

which is similar to Eq. (3.2). 
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Figure 4.2. a) Power spectrum density for the pressure data at buoy 1 without filtering, b) Power 
spectra density of the Hilbert transform after low pass filtering using a Chebyshev 
type II filter with a 0.08 cutoff frequency. 


Matlab’s functions angle and unwrap (phase in radians 2*pi wrap), provide 
the phase of the time series (4.1). Shown in Figure 4.3 are the graphs 
corresponding to the phase of the whole pressure data from buoy 1 and also for 
the two hours window. 
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Figure 4.3. Phase (rad) as a function of time (min) for the pressure data at buoy 1 for a) one day, 
b) two hours. 


The local frequency is by definition the slope of the phase and we applied 
to the phase data the MATLAB function diff and divided by the time between 
points (called tbp in Appendix) to get this slope. Figure 4.4 shows a plot of the 
local frequency as a function of time obtained from the data shown in Figure 4.3. 



Figure 4.4. Local angular frequency as a function of time for buoy 1 for a two hour time window. 
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Figure 4.5 shows the cross-correlation of data from buoys 1,2, and 3 
using the xcorrfunction in MATLAB. The cross-correlation is blown-up in Figure 
4.6. At the maximum of the cross-correlation we can obtain the time when the 


signals from two different buoys are highly correlated. Using this time, and 
knowing the distance between buoys, we can estimate the phase velocity at 
which the wave is moving from one buoy to the other. Assuming that the three 
buoys started to record the data at the same time, we get the phase velocities 
estimates shown in Table 4-1. 
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Figure 4.5. Cross correlation of the phase between a) buoys 1 and 2, b) buoys 1 and 3, and c) 
buoys 3 and 2. 
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Figure 4.6. Blown up of the cross correlations in Figure 4.5. 



distance 

correlation time 

calculated phase 
velocity 

Buoy 1 - Buoy 2 

123.46 m 

34 s 

3.63 m / s 

Buoy 2 - Buoy 3 

219.16 m 

163 s 

1.35 m / s 

Buoy 1 - Buoy 3 

342.62 m 

102 s 

3.36 m / s 


Table 4-1 . Relations of the distances among the buoys, correlation time and phase veiocity. 
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V. DISCUSSION OF RESULTS AND RECOMMENDATIONS FOR 

FUTURE WORK 


We have investigated the Hilbert transform as a tool to determine phase 
velocity from ocean wave data of three buoys. The values tabulated In Table 4-1 
were obtained by dividing the distance between buoys by the time for maximum 
phase correlation of the time series of pressure measurements at the 
corresponding buoys. 

Table 5-1 is a comparison between the values of phase velocity as 
determined by the correlation times to the values of phase velocity assuming 
linear shallow water theory over uniform bottom topography. The value of the 
phase velocity determined from shallow water theory represents an upper bound. 
If the normal to the wave front makes an angle with the line joining two buoys, 
the value of phase velocity obtained from correlations of the time series will yield 
an underestimation. Thus, the low values obtained from the correlations may be 
due, in part, to an effective average over all possible directions. 



phase velocity from 
correlations 

average depth 

calculated phase 
velocity 

Buoy 1 - Buoy 2 

3.63 m/s 

6.10 m 

7.73 m/sec 

Buoy 2 - Buoy 3 

1.35 m/s 

3.81 m 

6.11 m / sec 

Buoy 1 - Buoy 3 

3.36 m/s 

5.34 m 

7.23 m / sec 


Table 5.1 Comparison of phase velocity values determined from correlations (second 

column) and from shallow water theory (fourth column) assuming uniform depth 
(third column). 
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It is then apparent that to accurately determine the phase velocity from 
time series, we require an accurate directional spectral data or surface height 
time series from a bigger array of buoys. An accurate directional spectra 
provides the necessary angular information that can be used to determine the 
true value of the phase velocity. As shown by Rikiishi (1978) time series from a 
minimum of three buoys cannot be used to determine directional spectra for a 
broadband distribution of waves. 

Any data analysis on time series from pressure gauge measurernents 
have fundamental limitations. First, because the pressure gauge responds to 
the rise and fall to the free surface, the pressure at depth z attenuates by a factor 
proportional to e‘‘“ where k is 2%!% and X is the wavelength of the Fourier 
component of the surface disturbance. Unless the wave data corresponds 
exclusively to shallow water, the correlation time determined from the phase of 
the Hilbert transform at two points may be inaccurate. Second, if the wave field 
is nonlinear, the pressure and the surface height are nonlinearly related. Due to 
nonlinearities, the correlation time from pressure gauge measurements may yield 
misleading phase velocity results. Thus, phase velocity values from phase 
correlations using the Hilbert transform, are more reliable with time series of 
surface height measurements. 

If the deficiencies noted above can be corrected, the preliminary results of 
this thesis are evidence of the powerful technique of using Hilbert transform to 
determine the phase velocity of low frequency shoaling waves. This technique 
can be summarize operationally as follows. The analytic signal of a time series 
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possesses a natural phase. Correlations of the phase yield an average time of 
travel. If directional spectra are available, the angular information can then yield 
the true time of travel from which the phase velocity can be determined. In this 
approach, there is no explicit assumption of the relation between frequency and 
wavenumber of waves in the wave field. Indicating that it may be applicable to 
arbitrary wave fields, both linear and nonlinear. 


29 





30 




APPENDIX MATLAB CODE 


This appendix contains the particular mfile (Matlab Code) that we used to 
analyze the data from the three buoys in order to extract the information that we 
were interesting on. 

%Loading the desired data 
load 16jul1a.dat 
p1=16jul1a(:,1) 
clear 16jul1a 
load 16jul2a.dat 
p2=jul2a(:,1) 
clear 16jul2a 
load 16jul3a.dat 
p3=16jul3a(:,3) 
clear 16jul3a 
p1=p1-mean(p1) 
p2=p2-mean(p2) 
p3=p3-mean(p3) 
fs=2 

M1= max(size(p1)) 

N=172800 
f=0:fs/N:(N-1)*(fs/N) 
tbp=0.5 

t1=[1:M1].*(tbp/60) 

% Looking for the FFT of the data from the three buoys to select the 
desired frequencies 
P1=fft(p1)*2/length(p1) ; 

P2=fft(p2)*2/length(p2) ; 

P3=fft(p3)*2/length(p3) ; 

% Applying the filter and Hilbert Transform to the data 

cf_low=0.08 ; %cut off frequency for buoys 1,2 and 3. 

[b1,a1]=cheby2(5,40,cfJow/(fs/2));%Chebyshev filter 5 order 40 dB for data 

%form buoy 1 


%Sampling frequency 
%Get Number of data points 
%Total number of points 
%Normalizing the x-axis to frequency 
%Time between data points in seconds 
%Normalizing the x-axis to time in minutes 
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Y1 =filtfilt(b1, a1, p1) ; %Filtering data from buoy 1 

clear p1 ; %Clear data pressure 1 from memory 

V_wavep1 =hilbert(Y 1) ; %Hilbert Transform data from buoy 1 

phasepi =angle(V_wavep1) ; %Phase of the Hilbert Transform for buoy 1 

[S1 ,f 1 ]=psd(Y1,1024,fs) ;%examine spectrum of the low wave 

%without filtering 

[S2,f2]=psd(V_wavep1,1024,fs) ;%examlne spectrum of the low wave 

%after filtering and the Hilbert transform 
[b2,a2]=cheby2(5,40,cf_low/{fs/2)); %Chebyshev filter order 5,40 dB. 
Y2=filtfllt(b2, a2, p2) ;%Filtering data from buoy2 

clear p2 ; %Clear data pressure 2 from memory 

V_wavep2=hilbert(Y2) ; %Hilbert Transform data from buoy 2 

phasep2=angle(V_wavep2) ; %Phase of the Hilbert Transform for buoy 1 
[b3,a3]=cheby2(5,40,cf_low/(fs/2)); %Chebyshev filter for data from buoy 3 
Y3=filtfilt(b3, a3, p3) ; %Filtering data from buoy 3 

V_wavep3=hllbert(Y3) ; %Hilbert Transform data from buoy 3 

phasep3=angle(V_wavep3) ; %Phase of the Hilbert Transform for buoy 1 
clear p3 ; %Clear data pressure 3 from memory 

%Looking in a two hours windows of the data for each buoy 
V_wavep1 =V_wavep1 (1:14400) ; %points correspond to a two hours windows 
V_wavep2=V_wavep2(1:14400) ; %points correspond to a two hours windows 
V_wavep3=V_wave p3(1:14400) ; %points correspond, to a two hours windows 
Y1 =Y1 (1:14400) ; %points correspond to a two hours windows 

Y2=Y2(1:14400) ; %points correspond to a two hours windows 

Y3=Y3(1:14400) ; %points correspond to a two hours windows 

t1=t1 (1:14400) ; %Two hour time window 

%Looking for the phase of the Hilbert Transform for the two hour window 
theta1=unwrap(angle(V_wavep1)); %Phase of the Hilbert transform in radians 
theta2=unwrap(angle(V_wavep2)); %Phase of the Hilbert transform In radians 
theta3=unwrap(angle(V_wavep3)): %Phase of the Hilbert transform in radians 
%Looking for the derivative of the phase of the Hilbert Transform for the 
two hour window 

dthetal = diff(theta1 ./tbp) ; %By definition the Local frequency at 




cltheta2 = diff(theta2./tbp) 


%buoy 1 

;%By definition the Local frequency at 
%buoy 2 

dthetaS = diff(theta3./tbp) ;%By definition the Local frequency at 

%buoy 3 

%Looking for the cross-correlation of the phase of the Hilbert transform 
Cp21=xcorr(phasep2, phasep1,’coeff’) ;%Xcorrelations of the phase of the 

%slgnal from Buoy 2 and 1. 

Cp31=xcorr(phasep3, phasep1,’coeff’) ;%Xcorrelations of the phase of the 

%signal from Buoy 3 and 1. 

Cp32=xcorr(phasep3, phasep2,’coeff’) ;%Xcorrelations of the phase of the 

%signal from Buoy 3 and 2. 

Cps21 =fftshift(Cp21) ; %Shift of the correlations 

Cps31 =fftshift(Cp31) ; %Shift of the correlations 

Cps32=fftshift(Cp32) ; %Shift of the correlations 

%Plots of the FFT of the data pressure from the three buoys 
subplot (3,1,1) 

plot(f,abs(P1)),axis([0 0.5 0 0.05]),grld,title(’FFT of Data Pressure 
of buoy 1 vs Frequency’), 
subplot (3,1,2) 

plot(f,abs(P2)),axis([0 0.5 0 0.05]),grid, 
title(’FFT of Data Pressure of buoy 2 vs Frequency’), 
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subplot (3,1,3) 

plot(f,abs(P3)),axis([0 0.5 0 0.05]),grid, 

title(’FFT of Data Pressure of buoy 3 vs Frequency’), 

xlabel(’Frequency’),ylabel(’Normalized magnitude’) ; 

clear PI; clear P2; clear P3 ;%Release memory 

%Plot of the Power Estimate Spectrum 
figure; clf; 
subplot(2,1,1) 

loglog(f1,S1),axis([1E-3 1 IE-6 1]); 

title(’Power spectrum Estimate of Pressure data frequency without filtering’); 
ylabel(’magnitude (dB)’); 
subplot(2, 1, 2) 

loglog(f2,S2),axis([1E-3 1 IE-6 1]); 

title(’Power spectrum Estimate of V_wavep1 after Hilbert Transform ’); 
xlabel([’frequency(Hz)-cutoff frequency’,num2str(cfJow),’Hz’]); 
ylabel(’magnitude (dB)’); 

%Plot of the Phase in radians (Unwrap) of the Hilbert transform 

figure; clf; 
subplot(3,1,1) 
plot(t1,(phase_1)), 

title(’Plot of the phase of Hilbert Transform of the pressure data buoy 1’); 
xlabel(’time(min)’); ylabel(’Phase(rad)’); 


34 



subplot(3,1,2) 

plot(t1 ,(phase_2)),title(’Plot of the phase of Hilbert Transform of the pressure 
data buoy 2’); xlabel(’time(min)’);ylabel(’Phase(rad)’); 
subplot(3,1,3) 

plot(t1 ,(phase_3)),title(’Plot of the phase of Hilbert Transform of the pressure 
data buoy 3’); xlabel(’time(min)’); ylabel(’Phase(rad)’); 

%Plots of the Local frequency as a function of the time, data from buoy 1 

figure; clf; 

plot(t1,dtheta1); title(’Plot of the Local Frequency as a function of time at buoy 
1’); xlabel(’time(min)’); ylabel(’frequency(rad/sec)’); 

%Plots of the correlations of the data of the three buoys 

flgure;clf; 
subplot(3,1,1) 

plot(t1,abs(Cp21)),grid,title(’Correlation of Phase of the data from Buoy 2 and 
Buoy 1’),xlabel(Time(sec)’); 
subplot(3,1,2) 

plot(t1,abs(Cp31)),grid,title(’Correlation of Phase of the data from Buoy 3 and 
Buoy T),xlabel(Time(sec)’); 
subplot(3,1,3) 

plot(t1 ,abs(Cp32)),grid,title(’Correlation of Phase of the data from Buoy 3 and 
Buoy 2’),xlabel(Time(sec)’); 

%Plot of the shift of the correlations of the Phase from the three buoys 
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figure;clf; 
subplot(3,1,1) 

plot(t1,abs(Cps21)),grid,title(’Shift of Correlation of Phase of the data from Buoy 

2 and Buoy 1’),axis([0 1500 0 1]); 
subplot(3, 1,2) 

plot(t1 ,abs(Cps31)),grid,title(’Shift of Correlation of Phase of the data from Buoy 

3 and Buoy 1’),axis([0 1500 01]); 
subplot(3,1, 3) 

plot(t1,abs(Cps32)),grid,title(’Shift of Correlation of Phase of the data from Buoy 
3 and Buoy 2’),axis([0 1500 0 1]); 

%Plot of the Maximum values of the shifted correlations 

figure;clf; 
subplot(3,1,1) 

plot(abs(Cps21)),grid,title(’Max. value of the Correlation of the Phase of the data 
from Buoy 2 and Buoy 1 ’),axis([0 40 0.7 0.85]); 
subplot(3,1,2) 

plot(abs(Cps31)),grid,title(’Max. value of the Correlation of the Phase of the data 
from Buoy 3 and Buoy 1’),axis([75 150 0.5 0.7]); 
subplot(3,1,3) 

plot(abs(Cps32)),grid,title(’Max. value of the Correlation of the Phase of the data 
from Buoy 3 and Buoy 2’),axis([120 200 0.5 0.66]); 
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